Assessment of Precision and Accuracy of Brain White Matter Microstructure using Combined Diffusion MRI and Relaxometry

Joint modeling of diffusion and relaxation has seen growing interest due to its potential to provide complementary information about tissue microstructure. For brain white matter, we designed an optimal diffusion-relaxometry MRI protocol that samples multiple b-values, B-tensor shapes, and echo times (TE). This variable-TE protocol (27 min) has as subsets a fixed-TE protocol (15 min) and a 2-shell dMRI protocol (7 min), both characterizing diffusion only. We assessed the sensitivity, specificity and reproducibility of these protocols with synthetic experiments and in six healthy volunteers. Compared with the fixed-TE protocol, the variable-TE protocol enables estimation of free water fractions while also capturing compartmental T2 relaxation times. Jointly measuring diffusion and relaxation offers increased sensitivity and specificity to microstructure parameters in brain white matter with voxelwise coefficients of variation below 10%.


Introduction
Multi-compartment models are widely used in MRI studies of biological tissues.Compartmentalization allows for dividing a complex heterogeneous structure into simpler, homogeneous units.In the brain white matter (WM), there are at least four major compartments: intra-axonal space (IAS), myelin, extra-axonal space (EAS) and free water (FW).Each compartment is distinguished by unique diffusion profiles (Novikov et al., 2019) and T 2 relaxation times (Veraart et al., 2018).For instance, the T 2 relaxation time of myelin is 10-40 ms (Mackay et al., 1994), while T 2 of the white matter tissue is 70-100 ms (Stanisz et al., 2005), and T 2 of cerebrospinal fluid (CSF) is 500-2100 ms (Piechnik et al., 2009).Therefore, the signal of myelin water decays much faster with the echo time (TE) than the other compartments and is typically not observable for clinical diffusion MRI (dMRI) protocols.However, modeling either relaxation or diffusion properties independently presents challenges in parameter estimation, even for a simple twocompartment model (Jelescu et al., 2016a).Combining diffusion and relaxation holds tremendous promise, as these two modalities provide complementary information for every compartment (Veraart et al., 2018;Lampinen et al., 2023;Tax et al., 2021).
The FW compartment, though small in most WM voxels, gains significance at WM-CSF boundaries or in cases of edema.It is characterized by isotropic diffusivity and considerably longer T 2 relaxation times compared to the IAS or EAS.Quantifying FW fractions is of great interest for clinical applications, such as inflammation (Jelescu and Fieremans, 2023), schizophrenia (Pasternak et al., 2012) and Parkinson's disease (Planetta et al., 2016).Techniques have been proposed to model the FW compartment in addition to diffusion tensor imaging (DTI) (Pasternak et al., 2009), or to diffusion kurtosis imaging (DKI) (Collier et al., 2018).Yet, the inverse problem of estimating FW fractions solely from diffusion data is still ill-conditioned (Golub et al., 2021).Biophysical models of diffusion have routinely omitted FW to avoid these instabilities.To resolve this, multidimensional diffusion weightings have been demonstrated to facilitate solving FW fractions in a three-compartment system analytically (Reisert et al., 2019).By utilizing the differences in diffusion and T 2 relaxation times between FW and other compartments, the combined diffusion-relaxometry approach may enhance the sensitivity of multi-compartment models to FW, thereby improving the overall accuracy of parameter estimation (Anania et al., 2022).
The Standard Model (SM) of WM was proposed as an overarching framework (Novikov et al., 2019) that encompasses many previously proposed multi-compartment dMRI models (Kroenke et al., 2004;Jespersen et al., 2007;Assaf et al., 2004;Alexander et al., 2010;Fieremans et al., 2011;Zhang et al., 2012;Kaden et al., 2016;Reisert et al., 2017;Assaf et al., 2008;Jespersen et al., 2010a;Sotiropoulos et al., 2012;Reisert et al., 2017;Novikov et al., 2018).WM voxels are composed of a collection of fiber fascicles, which serve as the elementary building block for generating diffusion signals.Each fiber fascicle is assumed to consist of several non-exchanging compartments representing the IAS, EAS and FW.The overall dMRI signal is the weighted sum of the contributions from all bundles inside a voxel, each oriented according to the fiber orientation distribution function (fODF), see Fig. 1.The SM can be extended onto the relaxometry by assigning a T 2 relaxation time to each compartment, and varying the TE via a joint diffusionrelaxation protocol (Veraart et al., 2018;Tax et al., 2021;Lampinen et al., 2019).
An extended MRI protocol for diffusionrelaxation beyond conventional clinical protocols is necessary to overcome the degeneracy in parameter estimation.For instance, varying TE alters the relative weights between compartments with different T 2 values (Benjamini and Basser, 2016;Veraart et al., 2018;Ning et al., 2019;Lampinen et al., 2019).In addition, strong diffusion-weightings make dMRI measurements more sensitive to tissue microstructure (Novikov et al., 2019;Jespersen et al., 2010b).Multidimensional dMRI encodes diffusion along multiple directions simultaneously in the same diffusion weighting, thereby probing the response to a 3 × 3 B-tensor with rank B > 1 (Mitra, 1995;Lasič et al., 2014;Topgaard, 2017), and hence provides richer characterization of tissue microstructure than conventional dMRI encodings along a single direction (Coelho et al., 2019).
Our aim has been to optimize a comprehensive diffusion-relaxation MRI protocol within 30 minutes, combining multiple TE, high b-values and multidimensional diffusion encodings, and evaluate its sensitivity, specificity, and reproducibility.This work builds further upon the protocol op- timization by Coelho et al. (2022), via incorporating varying TE.Lampinen et al. (2020) previously designed a 15-minute acquisition protocol for diffusion-relaxometry by optimizing for the Cramér-Rao Lower Bound (CRLB) of the SM parameters.While the CRLB represents a theoretical lower bound of the variance for an unbiased estimator, most SM estimators are biased in conventional datasets.Machine learning (ML)-based estimators are not exempt from biases as they tend to push estimates towards the prior mean due to the limited obtainable information and SNR (Coelho et al., 2021), but they have shown to reduce overall estimation error significantly (Coelho et al., 2022).Thus, we adopted the root-mean-squarederror (RMSE) of an ML-based estimator (Reisert et al., 2017) as the objective function for our protocol design.
In the following sections, we will first explain the SM for multidimensional dMRI at varying TE.Next, we will present our framework of protocol optimization and the optimized variable-TE protocol (27 min).Then, we will compare the performance of the variable-TE protocol we proposed with its two subsets, a fixed-TE multidimensional diffusion protocol (15 min) and a 2-shell dMRI protocol (7 min) that is commonly used for DKI (Jensen et al., 2005).We will assess the sensitivity and specificity of each protocol using the recently proposed Sensitivity-Specificity Matrix (SSM) (Liao et al., 2023), and we will evaluate the reproducibility of each using the coefficient of variation (COV) from the scan and rescan of six healthy volunteers.Furthermore, we will assess the feasibility of estimating FW fractions using a 2-shell dMRI protocol along with b 0 measurements at multiple TEs (8 min).

Theory
2.1.T 2 relaxometry and multidimensional dMRI T 2 relaxometry has been widely used in conjunction with multi-compartment models to study volume fractions of myelin, axons, and FW (Mackay et al., 1994;Whittall et al., 1997).In a voxel with multiple water pools, each with its own T 2 relaxation, the overall signal fraction is: where s 0 is the proton density, each f i ≥ 0. In the continuum limit, Eq. ( 1) becomes a Laplace transform of the density spectrum of T 2 values.Thus, sampling MRI data at multiple TE can, in principle, enable the quantification of relaxation properties from different water pools.Although this framework has been widely applied, it is ill-conditioned, just like the inverse Laplace transform (Epstein and Schotland, 2008).Small fluctuations in the measured signals can lead to a wide range of reasonable outputs, making it unstable for the signal-to-noise ratio (SNR) of conventional datasets.dMRI provides a similar, yet complementary, tissue probe to relaxometry.Rather than generating signal decay due to the relaxation of the transverse magnetization, these measurements encode information about the random motion of water molecules.Conventional dMRI uses linear tensor encoding (LTE) to probe diffusion along a single direction within a given diffusion weighting.To probe tissue anisotropy, these measurements are performed along multiple directions by rotating the diffusion encodings.In contrast, multidimensional dMRI (Westin et al., 2016;Topgaard, 2017) simultaneously encodes diffusion along multiple directions within one diffusion weighting.This generalizes the scalar b-value to a B-tensor: (2) where g(t) is the effective diffusion gradient.The trace b = tr B is the conventional b-value.These additional experimental settings provide complementary information to LTE and their combination better informs biophysical modeling (Coelho et al., 2019;Reisert et al., 2019;Lampinen et al., 2019Lampinen et al., , 2023)).
In this work, we focus on axially symmetric Btensors, which can be represented as (Topgaard, 2017) where ĝ is the gradient direction (the symmetry axis) and β is the shape parameter.In particular, β = 1 for linear tensor encoding (one nonzero eigenvalue, LTE), β = − 1 2 for planar tensor encoding (two nonzero eigenvalues, PTE) and β = 0 for spherical tensor encoding (three identical nonzero eigenvalues, STE).In other words, the rank of the B-tensor determines the number of dimensions being probed by the diffusion gradient.
An axially symmetric diffusion tensor Building upon the fixed-TE protocol, the variable-TE protocol introduces shells with TE values of 62 ms, 78 ms, and 130 ms.A 3D plot visualizes these protocols within the acquisition parameter space, accompanied by a tree structure for clear representation of their relationship.In this plot, each dot is color-coded to represent a distinct shell.Notably, the 2-shell dMRI and fixed-TE protocol are distributed on a blue plane with TE = 92 ms.
lel and perpendicular diffusivities D ∥ and D ⊥ , and by a unit vector n defining its principal axis.When probing axially symmetric Gaussian compartments with axially symmetric B-tensors (3), the diffusion attenuation (assuming the unattenuated signal is normalized to 1) becomes where ξ = ĝ•n.

SM for diffusion-relaxometry
In the SM framework, WM voxels are modeled as a collection of fiber segments (fascicles), elementary bundles of aligned axons, see Fig. 1.The fiber fascicle contains three non-exchanging compartments: IAS, EAS and FW.Myelin is not visible at the TE used in clinical dMRI.
The IAS compartment models axons as zeroradius cylinders (sticks) with volume fraction f , relaxation time T a 2 and diffusivity D a along the axon direction.The EAS compartment represents hindered diffusion with relaxation time T e 2 and diffusivities parallel and perpendicular to axons D ∥ e and D ⊥ e , respectively.The FW compartment with fraction f w and relaxation time T w 2 has an isotropic diffusivity D w = 3 µm 2 /ms (water diffusivity at body temperature).
At clinically relevant experimental settings, tissue compartments behave approximately Gaussian.Thus, the signal response of a single fiber fasciclethe kernel -is a sum of signals from three Gaussian compartments weighted by their fractions and T 2 values.For an arbitrary diffusion-weighting and TE, the kernel is where f e = 1 − f − f w and ξ = ĝ•n is the scalar product between the symmetry axis of the kernel, n, and the symmetry axis of the B-tensor, ĝ.
In an imaging voxel, fiber segments are typically not aligned but are oriented according to a probability distribution P(n), the fODF.Therefore, the SM signal becomes the convolution of the kernel response and the fODF: where s 0 is the proton density and S 2 dnP(n) = 1.Eq. ( 6) has a general form of the spherical convolution that has been widely used in dMRI (Healy Jr et al., 1998;Tournier et al., 2004;Anderson, 2005;Tournier et al., 2007), where the assumptions on the physics of diffusion and relaxation determine the specific functional form (5) of the kernel.

Factorization in the spherical harmonics basis
The spherical harmonics (SH) basis is a natural choice for functions on a sphere.Here, we can write the fODF as: where Y ℓm (n) are the SH basis functions and p ℓm are the SH coefficients.Likewise, the SM signal can be written in the SH basis, where the convolution becomes a product: are the projections of the SM kernel onto the Legendre polynomials P ℓ (ξ) (nonzero m projections are zero due to the axial symmetry of the kernel).
To remove the dependence on the choice of physical basis in three-dimensional space, we define the rotational invariants of the signal and fODF (Novikov et al., 2019): (10) The fODF rotational invariants are normalized such that p 0 = 1 and 0 ≤ p ℓ ≤ 1.Here, p ℓ can become a metric of fODF anisotropy for any l > 0, and we typically choose p 2 because it has higher SNR than the other high-order terms.
Thus, we can represent all rotationally invariant information in the SM signal as: (11) Note that rotational invariants of odd ℓ are zero due to time-reversal symmetry of the Brownian motion.
Through this transform, we effectively compress the SM signal S(b, β, ĝ, TE) to a few rotational invariants without loss of information in terms of estimating the kernel parameters.

MRI protocol optimization
The MRI acquisition protocol determines the amount of microstructure information encoded into the MRI signals.Previous studies have suggested high b-values (Novikov et al., 2019;Jespersen et al., 2010b), multidimensional dMRI (Mitra, 1995;Lasič et al., 2014;Topgaard, 2017) and varying TE (Benjamini and Basser, 2016;Veraart et al., 2018;Ning et al., 2019;Lampinen et al., 2019) can provide additional information to conventional 2-shell dMRI protocols.Here, we design our protocol by incorporating all of these features.We optimize our protocol based on the RMSE of an ML-based estimator (see (Coelho et al., 2022) and Section 3.2 for more details) , where ϵ θi is the RMSE for a given SM parameter θ i , N train is the number of training samples for the ML-based estimator.Each RMSE term ϵ θi is normalized by the range of its parameter w i .
The full parameter set is defined as The T 2 value of the FW compartment is much longer than that of the IAS and EAS compartment and is set to 500 ms in our model.The objective function F is minimized by finding the optimal combinations of shells with (b, β, TE, N dirs ), standing for b-value, B-tensor shape, echo time and number of distinct gradient directions, respectively.The individual RMSE was computed analytically as a function of the shells, SNR and the moments of the training set distribution as in (Coelho et al., 2021).We minimized Eq. ( 12) using the Self Organizing Migrating Algorithm (Zelinka, 2016).
The total acquisition time was constrained within 27 minutes.The maximum b-value was set to 10 ms/µm 2 .Given the hardware limitations such as maximum gradient amplitude and slew rate, we calculated the minimum encoding time for a given {b, β} combination (Sjölund et al., 2015;Lampinen et al., 2020).Each shell was assumed to have diffusion encodings along gradient directions uniformly distributed on a sphere for optimal acquisition.Furthermore, we fixed conventional multishells (2-shell dMRI protocol) into our protocol together with the protocol used by Coelho et al. (2022) (fixed-TE protocol), see Fig. 2, to enable future comparisons.

Parameter estimation
Conventional maximum likelihood estimation (MLE) is unstable for SM parameter estimation at realistic SNR levels or with limited data.Reisert et al. (2017) first proposed the use of ML-based algorithms to estimate SM parameters, which have shown lower estimation error than conventional MLE.The ML-based estimator learns the mapping from synthetic MRI signals to model parameters from synthetic data generated from the model.Such mapping can be written as: where θ is the ML estimator, y i N i=1 are the noisy measurements (rotational invariants S l (b, β, TE) in our case), W is the degree of the polynomial regression, and a j1,j2,...,j N are the regression coefficients computed during training, which are a function of SNR.Coelho et al. (2021) demonstrated that a cubic polynomial regression is sufficient to approximate this mapping at realistic SNR levels and that more complex regressions, e.g.fully connected neural networks or random forests, lead to identical RMSE values.After generating signals based on Eqs.(5)(6), Gaussian noise was added to the signal, as we assumed the Rician floor was minimized by denoising the complex-valued data (Lemberskiy et al., 2019).During the inference, we selected the corresponding regression coefficients based on the estimated SNR (Veraart et al., 2016a).The code for the ML-based estimation is now publicly available in the Standard Model Imaging (SMI) toolbox (https://github.com/NYU-DiffusionMRI/SMI).
The performance of SMI, like other ML algorithms, is influenced by the nature of its training data.This impact is mediated by the acquired measurements and their SNR.In scenarios where comprehensive protocols are employed, the algorithm's performance is less susceptible to variations in training data.Conversely, in cases with limited protocols containing less information, the performance is more sensitive to changes in the training data.Here, we used a uniform distribution for training to minimize the impact of prior information, with lower bounds [0.05, 0.05, 1, 1, 0.1, 0, 50, 50] and upper bounds [0.99, 0.95, 3, 3, 1.2, 0.2, 150, 120] for the full parameter set e , f w , T a 2 , T e 2 }.In the context of two-shell dMRI protocols, the FW compartment is often omitted, as such protocols do not provide enough information to accurately retrieve parameters from all three compartments.To explore this limitation, we trained the SMI estimator under two conditions: excluding the FW compartment (2-shell woFW) and including it (2-shell wFW).In the 2- shell woFW scenario, the FW fraction was set to zero in all training samples.It is also important to note that fixed-TE and two-shell dMRI protocols are insensitive to T 2 relaxation times, as they are acquired for a single TE value.

Sensitivity-Specificity Matrix
Sensitivity and specificity are two key aspects for measuring the performance of model parameter estimation.In reality, the estimation of one parameter may be interfered by the others, creating spurious correlations and reducing specificity.To quantify sensitivity and spurious correlations of parameter estimation, we consider the Sensitivity-Specificity Matrix (SSM) in noise propagations whose elements quantify relative changes of an estimated parameter θj with respect to the actual change of parameter θ i .Here, angular brackets denote averaging over the distribution of ground truths (the test set).The normalization by the mean values is introduced for convenience, to make the off-diagonal elements dimensionless.Practically, we evaluated the SSM from a linear regression of the estimates θj with respect to ground truths θ i of all N θ parameters in a test set.The test set was sampled from a uniform distribution that covered the most probable range of the SM parameters θ with lower bound [0.3, 0.3, 1.5, 1.5, 0.4, 0, 50, 40] and upper bound [0.8, 0.8, 2.5, 2.5, 1, 0.15, 120, 100].In the simulation, the SNR was set at 40 for the b 0 measurements with a TE of 92 ms.
Ideally, the SSM should be an identity matrix.The deviation from an identity matrix shows the limitations of the parameter estimation.In an SSM, the diagonal elements quantify the sensitivity, while the off-diagonal elements quantify spurious correlations and the lack of specificity.The elements of SSM are typically between -1 and 1, similar to the correlation coefficient, which makes them easier to interpret than RMSE.

In vivo MRI data
After informed consent, six healthy volunteers, aged 24 to 60 years old (mean age: 33.8 ± 13.4 years), including 3 females and 3 males, underwent brain MRI on a 3T Siemens Prisma, using a 64-channel head coil for reception.The total scan duration was 27 minutes, details in Fig. 2. Due to the restricted availability of MRI sequences at the moment of scanning, we used the PTE shell with β of −0.5, while the optimization resulted in β value of 0.6.Maxwell-compensated asymmetric waveforms (Szczepankiewicz et al., 2019b) were employed in the acquisition protocols using an inhouse diffusion sequence that enables user defined gradient waveforms (Szczepankiewicz et al., 2019a).A non-diffusion-weighted image with reverse phaseencoding was acquired to correct for EPI induced distortions (Andersson et al., 2003).Imaging parameters: voxel size 2 × 2 × 2 mm 3 , bandwidth 1818 Hz/pixel, GRAPPA factor 2, partial Fourier 6/8.Subjects were scanned, repositioned and subsequently rescanned with the same MRI protocol.The b-values and directions were randomized to avoid drift bias (Szczepankiewicz et al., 2021).
Magnitude and phase data were reconstructed.This was followed by a process of phase estimation and unwrapping, which set the stage for denoising coil images (Lemberskiy et al., 2019;Veraart et al., 2016b).After denoising, the coil images were combined and then processed by the DESIGNER pipeline (Ades-Aron et al., 2018;Chen et al., 2023), including steps of denoising (Veraart et al., 2016b), Gibbs ringing correction (Lee et al., 2021), and correction for eddy current distortions and subject motion (Smith et al., 2004).The last step was performed separately for each unique combination of (β, TE).Finally, a rigid registration was computed to align all images to the subset with the largest number of DWI.
Regions of interest (ROI) were automatically segmented by nonlinear mapping of the fractional anisotropy (FA) map onto the WM label atlas of Johns Hopkins University (JHU) (Mori et al., 2005).FA maps for each subject were computed from a DKI fit based on the 2-shell subset of the data (Jensen et al., 2005).Mean values of major WM regions were extracted after excluding outliers that fell outside of the physical range of SM parameters.
The full variable-TE protocol was used to jointly estimate the full SM parameter set θ for diffusionrelaxometry.A subset of this protocol acquired at TE = 92ms and a further subset with b = 0, 1, 2 ms/µm 2 were used for estimating SM parameters excluding T 2 .The hierarchy of the variable-TE, fixed-TE and 2-shell dMRI protocol is illustrated in Fig. 2, where the latter is a subset of the former.When comparing the variable-TE results with the fixed-TE and 2-shell dMRI protocols, we weighed the fractions f and f w estimated from the variable-TE protocol with corresponding e −TE/T2 factors, and rescaled them to sum to 1.The rescaled fractions contained T 2 weighting and were denoted as f (TE) and f (TE) w for distinction.In addition to the above three protocols, we evaluated a clinically feasible 8-minute protocol by combining the 2-shell dMRI protocol (TE=92 ms) with b 0 of multiple TEs (62, 78, 130 ms).This protocol, featuring varying TE, was specifically chosen to assess the feasibility of accurately estimating FW fractions using a protocol that is practical for clinical applications.

Coefficients of variation
The scan and rescan of every subject were registered to one another for evaluating voxelwise reproducibility.Coefficients of variation (COV) were computed from the scan and rescan using the formula COV = ⟨σ/µ⟩, where the standard deviation was estimated as σ = √ π 2 |x 1 − x 2 |, and the mean µ was averaged over the scan and rescan of all available protocols.This ensures that the normalization for the same parameter across different protocols is consistent.The COV was first calculated for each voxel and then averaged over all the WM regions.We illustrated the COV for the six subjects in Fig. 6.The COV served as a reproducibility metric to assess the robustness of the acquisition protocol and parameter estimation.In addition, we employed the concordance correlation coefficient (CCC) to measure the consistency between two different protocols, and between the scan and rescan of the same protocol.

Protocol optimization
We incorporated several "orthogonal" acquisition tactics, such as high b-values, multidimensional dMRI and multiple TE, into our optimized protocol by exploring the optimal combination of (b, β, TE, N dirs ).The tree structure in Fig. 2 demonstrates that the fixed-TE protocol is a subset of the variable-TE protocol, and the 2-shell dMRI protocol is a subset of the fixed-TE protocol.In the optimized variable-TE protocol, the fixed-TE segment from Coelho et al. (2022) was fixed into our optimization, while the remainder of the protocol was obtained directly as a result of our optimization framework.
The 3D scatter plot of Fig. 2 illustrates the precise positioning of every shell in the acquisition space of (b, β, TE).The 2-shell dMRI protocol includes three LTE shells at 0, 1, 2 ms/µm 2 .Aside from that, the fixed-TE protocol consists of three shells: a high-b LTE shell at 6 ms/µm 2 , a low-b PTE shell at 2 ms/µm 2 and a low-b STE shell at 1.5 ms/µm 2 (Coelho et al., 2022).In the fixed-TE protocol, all six shells are acquired at a TE of 92 ms.The variable-TE protocol extends this by including three additional TEs: 62, 78, and 130 ms.Within these additional shells of the variable-TE protocol, all feature a low b-value (≤ 2 ms/µm 2 ), with only one being a PTE shell and the others LTE shells.

Sensitivity and specificity
Sensitivity and specificity are two key aspects of a biophysical model.Here, using the same ML parameter estimation algorithm and uniform prior distributions, we assessed the information content of the three protocols through their SSM (Fig. 3).There are slight differences in the spurious correlations (off-diagonal elements of SSM) between 2shell woFW and 2-shell wFW.Specifically, when the FW compartment is included in the training data, the correlation between the estimated D ⊥ e and the ground truth of f increases.Conversely, when FW is excluded, the estimated D ⊥ e shows a higher correlation with the ground truth of f w .
All parameters improve with additional "orthogonal" acquisitions, though this comes at the cost of increased acquisition time.For f and p 2 , which are very accurate already with 2-shell, the spurious correlations between them and compartmental diffusivities decrease by adding free gradient waveforms (fixed-TE) and get further reduced by adding multiple TE (variable-TE).D a and D ⊥ e estimations are substantially improved in fixed-TE and even more so in variable-TE.In particular, D a reaches a performance close to p 2 and f with variable-TE.On the contrary, D ∥ e shows the least improvement of accuracy after adding multidimensional dMRI and varying TE, indicating it is the least sensitive parameter.
The obvious gain from including varying TE is the ability to capture T a 2 and T e 2 , to which 2-shell and fixed-TE protocols are insensitive.The estimation of T a 2 is at the level of p 2 and f , and therefore very good, with the estimation of T e 2 showing a lower but still reasonable sensitivity.Yet, the most remarkable improvement happens to f w .In the 2-shell dMRI protocol, estimating the FW fraction is unfeasible.After including multidimensional dMRI, the estimated f w still exhibits strong spurious correlations with D ∥ e and D ⊥ e .In contrast, the variable-TE protocol demonstrates high sensitivity in estimating f w within a range of 0 to 0.15 in the test set, with the sensitivity metric at 0.99 (diagonal elements of SSM).The improvement for f w is expected since the signal contrast between FW and the remaining compartments increases when varying TE.

Parametric maps
The parametric maps of the three protocols and their voxelwise scatter plots are shown in Fig. 4. The parametric maps of f and D a show weaker contrast variation for the 2-shell protocol as compared to the other two richer protocols.In the scatter plots, we use the CCC (ρ) as the consistency metric to take into account the difference in the mean and standard deviation between two sets of estimations.Using estimates from variable-TE as the ground truth, as expected, fixed-TE has higher accuracy than 2-shell dMRI, and p 2 and f are overall more accurate than compartment diffusivities.Interestingly, 2-shell wFW has slightly higher consistency with variable-TE for f , D ∥ e and D ⊥ e than 2-shell woFW.Moreover, the disparity between variable-TE and its subsets highlights the extent of improvement gained through additional measurements.For instance, in the comparison between 2-shell wFW and variable-TE, D a (ρ = 0.06) and f (TE) w (ρ = 0.14) show the lowest correlations and thus the largest improvement after including additional data, followed by D ⊥ e (ρ = 0.30).When comparing fixed-TE with variable-TE, f (TE) w sees the largest improvement after adding shells with varying TE.Yet, the consistency levels between 2shell wFW, fixed-TE and variable-TE are relatively high for D ∥ e , suggesting minimal improvement for this parameter after acquiring more data.These findings are consistent with our previous analysis of the SSM.
Fig. 5 displays the ROI means of six major WM regions for the three protocols.Generally, the ROI mean differences between these protocols align with the voxelwise biases illustrated in the scatter plot of Fig. 4. The more advanced variable-TE protocol often yields estimates that deviate further from the mean of prior distributions, highlighting that the estimates are more informed by actual measurements rather than being dominated by the prior distribution.For instance, considering the prior mean for D a at 2 µm 2 /ms, variable-TE shows D a approximately 2.4 µm 2 /ms in the posterior limb of internal capsule (PLIC), while the 2-shell protocol yields D a close to the prior mean and fixed-TE shows intermediate values.Furthermore, the variable-TE protocol exhibits greater biological variability across different ROIs, while the 2-shell protocol may yield nearly uniform values for all ROIs, a trend that is particularly pronounced in the case of f w .In addition, using variable-TE as a reference, 2-shell wFW demonstrates smaller biases than 2-shell woFW for D ⊥ e and D ∥ e .

Reproducibility
Fig. 6 plots the voxelwise COV of each SM parameter based on a scan and rescan using the three protocols.Overall, the COV of most SM parameters are below 10%.The COV of 2-shell for f w is extremely low, and the variable-TE protocol significantly reduces the variation of f w between scan and rescan compared to fixed-TE.Fig. 7 presents the parametric maps of the scan and rescan using the variable-TE protocol.FA maps are plotted alongside SM parameters for reference.The CCC is calculated for all WM ROIs to evaluate the consistency between the scan and rescan.Parameters p 2 (ρ = 0.97) and f (ρ = 0.95) demonstrate the highest reproducibility, comparable to that of FA (ρ = 0.96).The reproducibility of all parameters in the variable-TE protocol are reasonably good, with the lowest CCC being above 0.7.

Free water estimation
Free water estimation is challenging with only 2shell protocols, and is still hampered by spurious correlations with EAS diffusivities using the fixed-TE protocol (Fig. 3).However, employing a variable-TE protocol considerably enhances the sensitivity to f w , prompting us to investigate the effectiveness of a clinically feasible protocol that combines a 2-shell dMRI protocol with b 0 measurements at three distinct TEs (62, 78, 130 ms) for FW estimation.Remarkably, this 8-minute protocol achieves sensitivity of 0.87 for f w and maintains spurious correlations below 0.30 (Fig. 8A).Although resolving other parameters like compartmental diffusivities and compartmental T 2 remains challenging, this 8-minute protocol, varying TE solely for b 0 measurements, is sensitive to FW fractions.Fig. 8B shows a scatter plot comparing f w estimates obtained from the streamlined 8-minute protocol against those from a more comprehensive 27-minute protocol.The CCC is 0.80 between these two protocols for f w , exceeding the CCC of 0.45 between fixed-TE and variable-TE shown in Fig. 4.This result highlights the benefit of varying TE for estimating free water fractions.

Discussion
In this paper, we presented an optimized 27-minute MRI acquisition protocol that extends conventional 2-shell protocols with three major features: high bvalue, multiple B-tensor shapes and multiple TE.This protocol enables a simultaneous measurement of both diffusion and relaxometry.We found the optimal protocol in the sampling space of (b, β, TE) by minimizing the RMSE of an ML-based estimator (SMI).This is different from the work from Lampinen et al. (2020) where CRLB was the objective function and MLE was adopted to estimate parameters.Our objective is tailored for the ML-based estimation, as the analytical solution of RMSE takes into account the distribution of training sets (Coelho et al., 2021).Then, we evaluated the sensitivity, specificity and reproducibility of the proposed 27-minute protocol and its subsets, a 7-minute 2-shell dMRI protocol and a 15-minute fixed-TE protocol, to demonstrate the gain in each model parameter from adding extra measurements.

Protocol optimization
In the optimization process, we integrated the fixed-TE protocol (Coelho et al., 2022) into our optimized design to enable comparisons.The final optimized protocol, which already includes multidimensional dMRI and high-b shells, was further enhanced by adding shells with variable TE.These added shells feature a TE range from 62 ms to 130 ms, compared to the uniform TE of 92 ms in the fixed-TE protocol.All of these additional shells have b-value between 0 and 2 ms/µm 2 , with 7 out of 8 additional shells being LTE.This configuration suggests that the variable-TE shells are primarily used to assess  signal decay with respect to TE in the low-b regime for LTE.

Sensitivity and specificity
Sensitivity and specificity of the different MRI protocols were quantified by the SSM.The SSM is conceptually similar to the confusion matrix in classification problems.In a classification problem, the prediction of a certain category may contain misclassified samples from other categories.Likewise, in a regression problem, the prediction of a certain parameter may contain information from other parameters.The spurious correlations undermine the core motivation behind biophysical modeling -specificity.Despite the effort to parameterize a quantity directly related to microstructure, estimating a nonlinear model parameter from limited noisy data can result in biases and spurious correlations.Constraints on the SM may yield relatively more precise results (Zhang et al., 2012;Kaden et al., 2016;Fieremans et al., 2011), but they do not resolve the issue of spurious correlations between model parameters (Liao et al., 2023).We demonstrate that acquiring "orthogonal" data is essential to address these issues effectively.
The 2-shell dMRI protocol of this study is largely the same as protocols adopted in some large imaging consortia, such as UK Biobank (Miller et al., 2016), Human Connectome Project (Glasser et al., 2016), Alzheimer's Disease Neuroimaging Initiative (Jack Jr et al., 2008) and Adolescent Brain Cognitive Development (Casey et al., 2018).Here, we study the gain of incorporating additional measurements by comparing 2-shell dMRI and fixed-TE with variable-TE protocols.With high b-value and multiple B-tensor shapes in the fixed-TE protocol, estimates of D a and D ⊥ e become much more sensitive and specific, while estimating f w and D ∥ e remain more challenging.Furthermore, varying TE improves the estimation of f w , which is anticipated due to the substantial difference between T 2 of FW and that of the IAS or EAS.The variable-TE protocol also reliably estimates T a 2 and T e 2 , both of which have the potential to become biomarkers for WM pathology.D ∥ e is improved in variable-TE, but caution is advised for its interpretation due to the relatively low sensitivity and specificity.
To estimate FW fractions for clinical applications, We designed a clinically feasible 8-minute protocol by combining the 2-shell dMRI protocol with b 0 measurements at three different TEs.This protocol does not require free gradient waveforms, and only extends the standard 2-shell protocol by adding several b 0 measurements of multiple TEs.This protocol, a small subset of the 27minute variable-TE protocol, demonstrates remarkably high sensitivity to f w (Fig. 8).It outperforms the 15-minute fixed-TE protocol for f w in terms of the SSM as well as the CCC when comparing against the 27-minute variable-TE protocol.Fig. 8B shows that the CCC is 0.80 for the entire range of f w , with even greater consistency observed for free water fractions over 5%, compared to the CCC of 0.45 for fixed-TE protocols against variable-TE for f w in Fig. 4.This indicates that varying TE can enhance f w estimation more effectively than using multidimensional dMRI.Nevertheless, this protocol falls short in accurately resolving compartmental diffusivities and T 2 values.When obtaining f w is of interest in clinical settings, a standard 2-shell dMRI protocol could be augmented with an additional minute of scanning time to include b 0 measurements at multiple TEs.

Reproducibility
We assessed the reproducibility of three protocols by scan/rescanning of six subjects.Overall, the COV of most parameters is between 5-10% voxelwise, which is comparable to FA.The COV of f w for fixed-TE protocols is high, in part due to the normalization by their small values as denominators.Yet, the reproducibility improved with variable-TE protocols, highlighting the benefit of incorporating multiple TE to estimate f w .Surprisingly, the 2shell dMRI protocol exhibits an extremely low COV for f w , as well as for some parameters that are typically difficult to estimate, such as D a .This somewhat unphysically low variation stems from the protocol's estimates being heavily influenced by the prior distribution.In light of that, the COV for ML-based estimators should be interpreted with caution and, ideally, in conjunction with the SSM.
The CCC is identified as an alternative metric for assessing reproducibility.In contrast to COV, which is calculated on a voxel-by-voxel basis before averaging over an ROI, CCC evaluates the entire dataset collectively.This characteristic makes CCC less susceptible to extreme values or small denominators in the normalization process.Besides, CCC accounts for the difference in the mean and standard deviation of two sets of measurements, going beyond just the Pearson correlation.As a standard metric for consistency, CCC is naturally suited for measuring reproducibility.Results from Fig. 7 suggest that p 2 and f are the only two parameters demonstrating a level of reproducibility comparable FA.This finding aligns with their estimation performance shown in the SSM (Fig. 3).

Training data
The selection of training sets, which serve as the prior distribution in this Bayesian approach, is crucial in regularizing ML-based estimators.These estimators are designed to minimize the overall RMSE across all training samples, rather than optimizing for each individual sample.In doing so, these estimators tend to strike a balance between bias and precision.For parameters that are more challenging to estimate, the estimator often shifts predictions toward the prior mean to enhance precision, albeit at the expense of increased bias.This approach aligns with Bayesian principles, where parameter predictions are based on a pre-existing belief about their distribution (Reisert et al., 2017).When estimating difficult parameters, adopting a more informative and narrow prior can enhance performance within the most probable parameter range (Liao et al., 2023).However, the inherent dependency of ML-based estimators on the prior distribution may introduce a systematic bias to the overall estimation if the SNR is too low or measurements are insufficient.Note that this effect is present in any data-driven estimator minimizing RMSE and is mediated by the SNR and the measurements acquired.Thus, to ensure consistency across different estimations, the same prior distribution needs to be used to train the ML-based estimator.
The bias due to the prior distribution is evident in the relatively flat parametric maps of the 2-shell dMRI protocol shown in Fig. 4. Specifically, for parameters like f w , which are challenging to estimate at low values using only the 2-shell protocol, the sensitivity and specificity are the lowest among all parameters.As a result, the ROI mean values of f w for the 2-shell protocol appear uniformly close to the prior mean.The lack of variation in f w also leads to an extremely low COV for 2-shell protocols between the scan and rescan.However, knowing that the FW fraction typically ranges between 5 and 10% from the parametric maps of variable-TE protocols, incorporating a similar amount of FW in the training samples for 2-shell protocols can mitigate bias in EAS compartment diffusivities when comparing 2-shell wFW with 2-shell woFW (Fig. 5).
As more "orthogonal" data is acquired, the dependence of estimates on the prior distribution decreases.ROI mean values from the variable-TE protocol tend to diverge more from the prior mean than those from less comprehensive protocols.While the 2-shell protocol shows both D a and D ∥ e to be around 2 µm 2 /ms (prior mean), the variable-TE protocol reveals that D a is consistently higher than D ∥ e across all examined ROIs.In addition, the variable-TE protocol captures a greater biological variability, as reflected in both the ROI means and the parametric maps.To enhance the sensitivity and specificity of SM parameter estimation, it is essential to incorporate additional measurements that can provide complementary information about the underlying tissue microstructure.

Parameter estimates in vivo
In addition to the simulation showing excellent performance (Fig. 3), the parameter estimates by the variable-TE protocol (Fig. 7) generally agree with the prior studies.The fraction of axons f shows the highest values in the most densely packed regions, such as corpus callosum and internal capsule.The mean value of IAS axial diffusivity D a is about 2.2 µm 2 /ms in WM, consistent with the result from planar signal filtering (Dhital et al., 2019).The FW fraction f w is very high near ventricles due to the partial volume effect of cerebrospinal fluid and fluid around the brain parenchyma.Otherwise, f w is between 0 and 5% in most WM voxels, and the T 2weighted FW fractions f (TE) w generally fall between 5 and 10%.T a 2 of the IAS is higher near ventricles also due to partial volume effect.T 2 ranges from 60 to 120 ms in the IAS and is generally higher than that of the EAS, which is consistent with findings in the human brain (McKinnon and Jensen, 2019;Veraart et al., 2018) and in ex vivo animal nervous tissues (Peled et al., 1999;Wachowicz and Snyder, 2002;Bonilla and Snyder, 2007;Dortch et al., 2010).

Limitations and outlook
In this section, we summarize the limitations and offer future outlook for this study.First, myelin water is characterized by a particularly short T 2 relaxation time (Mackay et al., 1994) compared to the other compartments of WM, as its signal decays rapidly after the MRI pulse.Due to the diffusion gradients employed in our MRI protocols, the TE of our measurements is beyond the range for capturing myelin water fractions.However, we expect the presence/absence of myelin to impact D ⊥ e (Jelescu et al., 2016b).D ⊥ e can serve as a proxy for myelin water fractions.Second, when comparing the three protocols, different models were employed.There are no T 2 relaxation parameters in the 2-shell dMRI and fixed-TE protocol, as they are acquired at only one TE.However, after reweighing fractions with their corresponding compartmental T 2 , f (TE) and f (TE) w from the variable-TE protocol achieve excellent agreement with those estimated from fixed-TE (Fig. 4).This result provides an important consistency check for SM assumptions, as the fractions and T 2 relaxation times can maintain their relationship when estimated separately.Third, we did not explore the optimization landscape or the protocol parameter in detail in this manuscript.Readers can refer to Coelho et al. (2022) for more discussions.The optimization landscape as a function of protocol parameters is quite flat around this protocol, and that means small changes in the shell settings will hardly affect the ultimate performance.Fourth, the time dependence of diffusion is ignored in this study, as we assume that for our experimental settings diffusion in WM has reached the longtime limit and thus every compartment can be modeled as Gaussian.Fifth, even in the variable-TE protocol, D ∥ e still has relatively low sensitivity and needs to be improved further for obtaining a robust estimation of the EAS axial diffusivity.Sixth, the 27-minute variable-TE protocol is time-efficient for estimating the entire set of SM parameters, as the optimization framework was aimed for.However, to obtain a specific parameter, a less comprehensive protocol that is practical for clinical settings might be adequate.The 8-minute protocol extending 2-shell protocols with b 0 of multiple TEs can achieve fairly high sensitivity to f w .Future studies can focus on optimizing the protocol for a selection of SM parameters, thereby addressing specific clinical needs.

Conclusion
In this work, we presented a protocol optimization framework that generated a 27-minute variable-TE protocol for joint modeling of diffusion-relaxation.The variable-TE protocol was compared to its subsets, a fixed-TE protocol and a 2-shell dMRI protocol.The results show that by varying TE, FW fractions can now be accurately estimated even at very low values, which is challenging with the fixed-TE protocol.In addition, compartment T 2 values are now simultaneously estimated with high accuracy.In terms of reproducibility for variable-TE protocols, COV are below 10% for almost all SM parameters, comparable to FA.Our work demonstrates the benefits of multidimensional dMRI and varying TE to the parameter estimation.The variable-TE protocol can become a powerful tool for studying brain tissue microstructure.

Figure 1 :
Figure 1: Diagram of the Standard Model of diffusion-relaxation in white matter.In an MRI voxel, multiple elementary fiber fascicles coexist (whose MRI signal yields the so-called fiber response or kernel).Fiber fascicles contain multiple compartments that are characterized by their diffusivities, relaxivities, and water fractions.An MRI voxel is thus modeled as a collection of such fascicles oriented following an arbitrary fODF P(n).

Figure 2 :
Figure 2: Hierarchy of the optimized MRI acquisition protocols.Each distinct shell is represented in the form of (b, β, TE, N dirs ), where b denotes b-value, β B-tensor shape, TE echo time, N dirs number of gradient directions.The broader fixed-TE dMRI protocol encompasses the 2-shell dMRI protocol, as well as a low-b PTE and a STE shell, in addition to a high-b LTE shell.Building upon the fixed-TE protocol, the variable-TE protocol introduces shells with TE values of 62 ms, 78 ms, and 130 ms.A 3D plot visualizes these protocols within the acquisition parameter space, accompanied by a tree structure for clear representation of their relationship.In this plot, each dot is color-coded to represent a distinct shell.Notably, the 2-shell dMRI and fixed-TE protocol are distributed on a blue plane with TE = 92 ms.

Figure 3 :
Figure 3: Sensitivity-specificity matrix of the 3 protocols.Diagonal elements measures the sensitivity of the estimation, while off-diagonal elements quantify (lack of) specificity.According to the SSM, the fixed-TE protocol significantly improves the estimation of Da and D ⊥ e , and the variable-TE protocol can further accurately estimate FW fractions as well as capture compartment T 2 values.

Figure 4 :
Figure 4: Parametric maps of SM parameters for the 3 protocols.Parametric maps of a 32-year-old male are shown for each protocol overlaid on the b 0 images.The mask for parametric maps is FA > 0.2.Scatter plots of voxelwise values from all WM ROIs are shown below with the concordance correlation coefficient ρ. f (TE) and f (TE) w are reweighed by estimated T 2 to compare with the T 2 -weighted fractions in the 2-shell dMRI and fixed-TE protocols.

Figure 5 :
Figure 5: ROI means of SM parameters for the 3 protocols.Mean values of major WM regions from 6 subjects are shown for different protocols.T 2 relaxation times are shown for the variable TE protocol alone, and FA values are calculated using only the 2-Shell protocol.Here, fractions f (TE) and f (TE) w are T 2 -weighted, hence the notation TE in superscript.The error bar of the bar plots indicates the standard deviation of mean ROI values among the 6 subjects.ROI abbreviations: CC, corpus callosum; ALIC, anterior limb of internal capsule; PLIC, posterior limb of internal capsule; ACR, anterior corona radiata; SCR, superior corona radiata; PCR, posterior corona radiata.

Figure 6 :
Figure6: Reproducibility of SM parameter estimation.The coefficient of variation (COV) is mostly below 10% for almost all SM parameter.The COV of FW fractions for 2-shell is extremely low due to the dominance of the prior distribution in parameter estimation, while the COV of fixed-TE is notably high because of the small values of FW fractions in the normalization of COV.Note that fractions for the 2-shell dMRI and fixed-TE protocol are T 2 -weighted, and fractions for the variable-TE protocol are not.We use the same symbol for every protocol for clarity.

Figure 7 :
Figure 7: Parametric maps of SM parameters for the variable-TE protocol.Parametric maps of a 28 year-old male are shown for the scan and rescan of variable-TE protocol overlaid on b 0 images.The maps are masked by FA > 0.2.Scatter plots of voxelwise values from all WM ROIs are shown below with the concordance correlation coefficient ρ.

Figure 8 :
Figure 8: Free water estimation using the 2-shell protocol augmented with b 0 of varying TE.Panel A shows the SSM for a clinically-feasible 8-minute MRI protocol, which integrates the 2-shell dMRI protocol (TE=92 ms) with b 0 measurements at varying TE (62, 78, 130 ms), specifically designed to estimate FW fractions.Panel B presents a scatter plot comparing fw estimated from this 8-minute protocol against those from a more extensive 27-minute variable-TE protocol, utilizing data from the same 32-year-old subject previously shown in Fig. 4. The concordance correlation coefficient ρ between these two protocols is indicated on the scatter plot.Please note that, here, fw is not T 2 -weighted, different from f (TE) w shown in Fig. 4. Panel C displays fw parametric maps of this subject generated by both protocols, illustrating the spatial distribution of FW content across the WM.